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Max-stable processes for modelling extremes observed in space 

and time 

Richard A. Davis * Claudia Kliippelberg ^ Christina Steinkohl ^ 

Abstract 



Max-stable processes have proved to be useful for the statistical modelling of spatial 
f\j ' extremes. Several representations of max-stable random fields have been proposed in the 

^^ ' literature. For statistical inference it is often assumed that there is no temporal dependence, 

i.e., the observations at spatial locations are independent in time. We use two representations 
M-l , of stationary max-stable spatial random fields and extend the concepts to the space-time 

domain. In a first approach, we extend the idea of constructing max-stable random fields 
as limits of normalized and rescaled pointwise maxima of independent Gaussian random 
fields, which was introduced by Kabluchko, Schlather and de Haan |23j . who construct max- 
stable random fields associated to a class of variograms. We use a similar approach based on 
a well-known result by Hiisler and Reiss [H] and apply specific spatio-temporal covariance 
^ I models for the underlying Gaussian random field, which satisfy weak regularity assumptions. 

' Furthermore, we extend Smith's storm profile model |38| to a space-time setting and provide 



^^ ■ explicit expressions for the bivariate distribution functions. 

The tail dependence coefficient is an important measure of extremal dependence. We 
show how the spatio-temporal covariance function underlying the Gaussian random field can 
be interpreted in terms of the tail dependence coefficient. Within this context, we examine 
different concepts for constructing spatio-temporal covariance models and analyse several 
specific examples, including Gneiting's class of nonseparable stationary covariance functions 
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1 Introduction 

The statistical modelling of extremes is an important topic in many applications. Environmental 
catastrophes like hurricanes, floods and earthquakes can cause substantial damage to structures 
like bridges, towers and buildings. Within an insurance context, those extremal events can result 
in large losses and it is essential to have an estimate of the risk of such catastrophes. Adequate 
stochastic models for the characterization and quantification of the behaviour of extremal events 
are needed, leading us to extreme value theory. 

Univariate extreme value theory is well developed and standard introductions can be found 
in many books, including for example Embrechts, Kliippelberg and Mikosch [IT], Coles [6] and 
Leadbetter, Lindgren and Rootzen [25]. The extension of the univariate to the multivariate 
extreme value distribution is considered for instance in de Haan and Resnick [15] and Beirlant 
et al. [3]. In the simplest case, the process underlying the extremes is assumed to be independent 
or at least stationary and satisfying a mixing condition (see Leadbetter [M|)- The development 
of nonstationary models for extremes is still in an evolutionary stage. In the literature covariates 
are often introduced to overcome the nonstationarity of the underlying process. See for instance 
Davison and Smith [12] and Coles [6], Chapter 6, for more details. 

Max-stable processes are a natural extension of multivariate extreme value distributions to 
infinite dimensions. Detailed introductions and different spectral representations of stationary 
max-stable processes have been developed for example in Deheuvels [16], de Haan [13], de Haan 
and Pickands [13] and Schlather |34j . So far, max-stable processes have mostly been used for the 
statistical modelling of spatial data. Examples are given in Coles [2] and Coles and Tawn [7], 
who model extremal rainfall fields using max-stable processes. Another application to rainfall 
data can be found in Padoan, Ribatet and Sisson [30], who also describe a practicable pairwise 
likelihood estimation procedure. An interesting application to wind gusts is shown in Coles and 
Walshaw [8], who use max-stable processes to model the angular dependence for wind speed 
directions. 

Different approaches for constructing max-stable processes have been introduced, and we 
mention two examples here, which we will extend to the space-time setting. In some cases the 
approaches can lead to the same finite-dimensional distribution functions. 

The idea of constructing max-stable random fields as limits of normalized and scaled point- 
wise maxima of Gaussian random fields was introduced in Kabluchko, Schlather and de Haan [23], 
who construct max-stable random fields associated with a class of variograms. The limit field in 
this approach has the same finite dimensional distribution as the above described Brown- Resnick 
process (see Brown and Resnick [3] or Kabluchko et al. [23]). 

In an earlier paper Smith [38] introduced another max-stable process, which became known 
as the storm profile model. The different variables in the construction have an interpretation as 
components in a storm, including the shape and the intensity. The process is based on points 
from a Poisson random measure {{^j,Uj),j = l,...} together with a kernel function /, which 
in particular can be a centered Gaussian density. The max-stable process then arises from 



rj{y) = maxj>i ^jf{y; Uj). 

In real world applications, measurements are typically taken at various locations, sometimes 
on a grid, and at regularly spaced time intervals. For statistical inference, it is then often as- 
sumed, that the measurements are independent in time. We mention two approaches proposed in 
the literature concerning the analysis and quantification of the extremal behaviour of processes 
observed both in space and time. A first idea for modelling extremes in space and time can be 
found in Davis and Mikosch [lOj . who study the extremal properties of a moving average process, 
where the coefficients and the white-noise process depend on space and time. Sang and Gelfand 
[33] propose a hierarchical modelling procedure, where on a latent stage spatio-temporal depen- 
dence is included via the parameters of the generalized extreme value distribution. Extremes of 
space-time Gaussian processes have been studied in Kabluchko [22]. He analyses processes of 
the form supj/grg^t^i Z{snS,t') for some suitable chosen space-time Gaussian process and shows 
that the finite dimensional distributions of a properly scaled version converge to those of a 
Brown- Resnick process. 

For the extension of Kabluchko et al.'s approach to a space-time model we need an underlying 
spatio-temporal correlation model for the Gaussian random field, which satisfies a certain regu- 
larity condition at 0. This condition, taken from a well-known result by Hiisler and Reiss [21] . 
assumes that a rescaled version of the correlation function p{h, u) for the space-lag h and the 
temporal lag u satisfies logn(l — p{snh,tnu)) — )• 6{h,u), as n — ;• oo for the scaling sequences s„ 
and tn- We establish an explicit connection between the limit 6{h,u) and the tail dependence 
coefficient for two locations at two time points. 

Recently, the development of covariance models in space and time has received much at- 
tention and there is now a large literature available for the construction of a wide-range of 
spatio-temporal covariance functions. Examples can be found in Gneiting [20] . Ma [281 126] . 
Cressie and Huang [9] and Schlather [35]. Within this context, we generalize an assumption on 
correlation functions for the analysis of extremes from stationary Gaussian processes (see for 
instance Leadbetter et al.[25]. Chapter 12), which is a sufficient condition for the limit assump- 
tion described above. We show how Gneiting's class of nonseparable and isotropic covariance 
functions [20] fits into this framework. 

In addition, we examine spatial anisotropic correlation functions, which allow for directional 
dependence in the spatial components. Perhaps the easiest way to introduce anisotropy in a 
model is to use geometric anisotropy. For a detailed introduction, we refer to Wackernagel [39] . 
Chapter 9. Using this concept in the underlying correlation function, we can model anisotropy 
in the corresponding max-stable random field. Furthermore, we revisit a more elaborate way of 
constructing anisotropic correlation models based on Bernstein functions introduced by Porcu, 
Gregori and Mateu [52]; called the Bernstein class. 

Our paper is organized as follows. Two approaches for the construction of space-time max- 
stable random fields are described in Section 2. The max-stable random field, which is based on 
the maximum of rescaled and transformed replications of Gaussian space-time random fields is 
developed in Section [2.11 while Smith's storm profile model is extended to a space-time setting 



in Section 12.21 In Section [3l we show how Pickands dependence function and the tail depen- 
dence coefficient relate to the correlation model used in the underlying Gaussian random field. 
Further correlation models are discussed in Section [J] and simulations based on a set of different 
parameters are visualized. Section [4.21 analyzes anisotropic correlation functions, where one can 
see directional movements in the storm profile model, which are not possible in the isotropic 
case. 

2 Extension of extreme spatial fields to the space-time setting 

Max-stable processes form the natural extension of multivariate extreme value distributions to 
infinite dimensions. In the literature, different approaches for establishing space-time max-stable 
processes have been considered. In this section, we discuss two approaches for constructing such 
random fields and extend the concepts to a space-time setting. In Section 12.11 we describe the 
construction introduced in Kabluchko, Schlather and de Haan [23] based on a limit of pointwise 
maxima from an array of independent Gaussian random fields. Furthermore, we extend the 
approach introduced in de Haan [13] and interpreted by Smith [38] as the storm profile model 
to a space-time setting in Section 12. 2[ 

2.1 Mcix-stable random fields based on spatio-temporal correlation functions 

Before presenting the construction of a max-stable Gaussian random field in space and time, we 
begin with the definition of the Brown-Resick space-time process with Frechet marginals (see 
Brown and Resnick 0] or Schalther [33]). Let {S,j,j > 1} denote points of a Poisson random 
measure on [0,oo) with intensity measure E,~'^d^ and let Yj{s,t),j = 1,2,... be independent 
replications of some space-time random field {Y{s,t),{s,t) G M'^ x [0, oo)} with E(y(s,i)) < oo, 
and Y{s,t) > a.s., which are also independent of S^j. The random field, defined by 

oo 

r,{s,t) = \/ ^,Yj{s,t), (s,i)GM'^x [0,oo) (2.1) 

i=i 
is a max-stable random field with Frechet marginals and often refered to as the Brown-Resnick 
process (see Kabluchko et al. [23]) The finite-dimensional distributions can be calculated using a 
point process argument as done in de Haan [13]. For example, if (si, ti), . . . , (sa', t/c) are distinct 
space-time locations (duplicates in the space or the time components are allowed), then 

Pivisi, h) < yi, . . . , v{sK, tx) <yK) = P (ij V ^^^^^ < 1, Vi = 1, 2, . . . I 

= P {N{A) = 0) = exp I -E (" y ^^^^) } ' (2.2) 

where A = {{u,v)]uv < 1} and A^ is the Poisson random measure with points at 

Yj{sk,tk ) 
Vk 




As we will see below, the Brown-Resnick process will be the limit of a sequence of pointwise 
maxima of independent Gaussian space-time processes. In the following, let {Z{s,t)} denote a 
space-time Gaussian process on M.'^ x [0, oo) with covariance function given by 

Cisuh; S2M) = Cov iZ{si,ti),Z{s2,t2)) , 

for two locations si,S2 G M'^ and time points ti,t2 G [0,oo). We assume stationarity in space 
and time, so that we can write 

C{si,ti; S2,t2) = C{si - S2,ti - t2) = C{h, u), 

where h = si — S2 and u = ti — t2- Furthermore, let p{h,u) = C{h,u)/C{0,0) denote the 
corresponding correlation function. We will assume smoothness conditions on p{-, •) near (0,0). 
This assumption is natural in the context of spatio-temporal random fields, since it basically 
relates to the smoothness of the underlying random field in space and time. 

Assumption 2.1. There exists two nonnegative sequences of constants s^ — ?• 0, i^ ~^ as 
n — 7- oo and a nonnegative function 6 such that 

logn(l - p{sn{si - S2),tn{ti - ^2))) "^ S{si - S2,*i " ^2) G (0,oo), n -^ OO, 

foraU{si,ti) /(S2,t2), (si,ii), (s2,i2) G M'>' x [0,oo). 

Examples of such correlation functions are given in Section^ If (si,ti) = (52,^2)5 it follows 
that the correlation function equals one, p(0, 0) = 1, which directly leads to (^(0,0) = 0. The 
following theorem regarding limits of finite-dimensional distributions stems from Theorem 1 in 
Hiisler and Reiss [21] . In the following let C{M.'^ x [0, 00)) denote the space of continuous functions 
on M.'^ X [0,00), where convergence is defined as uniform convergence on compact subsets K of 
M^ X [0,00). 

Theorem 2.2. Let Zj(s,t), j = 1,2, ... be independent replications from a stationary space- 
time Gaussian random field with mean 0, variance 1 and correlation model p satisfiying Assump- 
tion \2.1\ with limit function 6. Assume there exists a metric D onW^ x [0, 00) such that 

6isi - S2,t, - t2) < {D{{Si,ti), {S2, t2))f (2.3) 



and set 

rjr,{s,t) = -\J - (s,t)GM'^x[0,oo). (2.4) 

n .V^ log {9{Zj{SnS,tnt))) 

Then, 

r]n{s,t) — >ri{s,t), (2.5) 

where — > denotes weak convergence in C{W^ x [0, 00)) and |r/(s, t), (s,i) G M'^ x [0,oo)} is a 
space-time max-stable process. The bivariate distribution functions for ??(s,i) have an explicit 
form given by 

r 1 / log ^2 \ 1 / log 21 \ ^ 

Fiyuy2) = exp ci> -^=m + Vm;uj - -^ w== + V^(J^) ■ (2-6) 

I yi \2^d{h,u) I y2 \2^6{h,u) \ 



Remark 2.3. Condition ()2.ip is sufficient to prove tightness of the sequence {r]n)neN i^ C{M. x 
[0, oo)). 74s shown in the proof of Theorem 17 in Kahluchko et al. J23^ the limit process tj turns 
out to be a Brown- Resnick process with Y in (j2.ip given by 

exp {W{s, t) - 6(3, t)} , (s, t) G M'^ X [0, oo), 

where {W{s,t), (s,t) G M*^ x [0, oo)} is a Gaussian random field with mean and covariance 
function 

Cov{Wisi,ti), Wis2,t2)) = S{si,tl) + 6{s2,t2) - 6{SI - S2,ti - ts). (2.7) 

In particular, 6 is a variogram leading to a valid covariance function in (j2.7p . 

Proof. Although this proof is similar to the one given in Kabluchko et al. [23], we provide a 
sketch of the argument for completeness. We start with the bivariate distributions. Prom classical 
extreme value theory (see for example Embrechts, Kliippelberg and Mikosch ^^, Example 

3.3.29), we have for 

, /7T1 loglogn + log(47r) 

bn = V 2 log n (2.8) 

iyilogn 

that 

lim $" L + i^) = e-yy. 

By using the standard arguments as in Embrechts et al. [T7], it follows that 

By applying this relation and Theorem 1 in Hiisler and Reiss [2Tj to the random variables 
r]n{si,ti) and r/„(s2,i2) for fixed (si, ti), (^2,^2) G M*^ x [0,oo), we obtain 

P{r]nisi,ti) < yi,r]n{s2,t2) < y2) 

/ n J "1 '^ 

= ^ \y^~log{Z,{SnSl,tnt,)) - ''^''V"log(Z,(s2,tnt2)) " ""^^ 

(n n 

\l Z,{SnSl,tntl) < <P~' (e-l/("^i)) , V Z,{s2,tnt2) < $"^ (e^l/^"^^)) 

~ P ( Zi(s„Si,t„ti) < Vbn,Zi{SnS2,tnt2) < 7 ^ "n 

\ bn On / 

/ 1 1| uf^f ^ -^ \ ^ fog(2/i) , , 7 / + + ^ ^ fog(y2) , , 

~ exp<^ \-nP [ Zi(s„Si,t„ti) > \-bn,Zi{SnS2,tnt2) > 7 ^ »" 

[ yi y2 \ bn bn 



The vector (Zi(s.„si,i„ti), Z2(s„S2, tni2)) is bivariate normally distributed with mean and 
covariance matrix given by p(s,„(si — S2),tn{ti — ^2)). Using the properties of the conditional 
normal distribution and Assumption I2.H it can be shown that the last expression converges 



to ()2.6p . Similarly to the procedure above, the finite-dimensional limit distributions of beyond 
second order can be calculated by using Theorem 2 in Hiisler and Reiss |21j . 

It remains to show that the sequence {r]n) is tight in C(M"' x [0,oo)). Following Kabluchko 
et al. [23], the main step of the proof is to show that the conditional family of processes 
{y„"(s,t), is,t) G M'^ X [0,oo)}, given by 

Y^{s,t) = {bn{Z{SnS,tnt)-bn)-io)\{bn{Z{O,0)-bn = io), WG[-C,c], UGN 

is tight in C{K), where K is any compact subset of M x [0, oo) and 6„, is defined in (j2.8p . This 
is achieved by calculating an upper bound for the variance of the distance between the process 
at two spatio-temporal locations, which in our case is given by Assumption (j2.3p . That is, for 
large n, 

Var(y„"(si,ti) -y„"(s2,t2)) < 2blil-p{snisi - S2) , Uh - t2))) 

< 2B6{si - S2, h - t2) < 2BD{{si,h), {si,t2)f, 

where i? > is some constant. The rest of the proof follows analogously to the proof in [23]. D 

Remark 2.4. Kabluchko f2^ studies the limit behavior of rescaled space-time processes of the 
form 

sup Z{SnS,t'), 
t'e[0,in] 

and shows that a rescaled version converges in the sense of finite- dimensional distributions to a 
space-time Brown-Resnick process. The assumptions on the covariance function in the underlying 
Gaussian space-time random field are similar to those we use in Section^ The approach differs 
from ours in the sense that we analyse the pointwise maxima of independent replications of 
space-time random fields, rather than the supremum over time of a single random field. 

Remark 2.5. In applications, the marginal distributions are often fitted by a generalized extreme 
value distribution and are then transformed to standard Frechet. Sometimes it may be useful to 
think about other marginal distributions, such as the Gumbel or Weibull. In order to use Gumbel 
marginals, we need 

n 

r]nis,t) = y -log{-log{<P {Zj{snS,tnt)))) -login), (2.9) 

and obtain the bivariate distribution function in (|2.6|) with 1/yi and l/y2 replaced by e~^^ and 
g-j/2 ^ jj yjQ want to use Weibull marginals in our model, we obtain 

n 
r]n{s,t) =n\J log($ {Zj{SnS,tnt))) , (2.10) 

i=i 

leading to the same bivariate distribution function as in (j2.6p . but with 1/yi and 1/2/2 replaced 
by yi and 2/2. respectively. 



2.2 Extension of the storm profile model 

In this section, we extend the following max-stable process, first introduced in de Haan [13], to 
the space-time setting. The process was interpreted by Smith [38] as a model for storms, where 
each component can be interpreted as elements of a storm, like intensity or center. In later 
papers, including for instance Schlather and Tawn [36] this process is called the storm profile 
model. We extend the concept to a space-time setting, where extremes are observed at certain 
locations through time. For simplicity of presentation we assume without loss of generality that 
M? is the space domain. Assume, that we have a domain for point processes of storm centres 
Z CM."^ and a time domain X C [0, oo), for which the storm is strongest at its centres. Further, 
let {{S,j,Zj,Xj),j > 1} denote the points of a Poisson random measure on (0, oo) x Z x X with 
intensity measure ^~^d^xA2((i2)x Ai (fix), where A^ denotes Lebesgue measure on R for d = 1,2. 
Each ^j represents the intensity of storm j. Moreover, let f{s, t; z,x) for [s, t) G M x [0, oo) and 
{z,x) £ Z X X he a non- negative function with 



f{s,t;z,x)X2idz)Xi{dx) = 1, {s,t) G E^ x [0,oo). 

ZxX 

The function / represents the shape of the storm. Define 

V{s,t)=\/ {Cjfiz,,Xj;s,t)}, (s,t)GM2x [0,oo). (2.11) 

The product S,jf{zj,Xj; s,t) can be interpreted as the wind speed at location s and time point 
t from storm j with intensity ^j, spatial location of the center Zj and maximum wind speed at 
time Xj at the centre. The finite dimensional distribution function of (?7(si,ti), . . . ,ri{sK,tK)), 
defined for fixed (si, ii), . . . , {sk, tx) G M'^ x [0, oo) and yi, . . . ,yK G R, is given through the 
spectral representation calculated in de Haan [13] as 



F{y,,...,yK) = expi- I V^^^^^^f^^^2(d^)Ai(dx)l. (2.12) 

To connect the storm model with the Brown- Resnick process that arises in Theorem 12.21 we 
assume a trivariate Gaussian density for the function / with mean (z, x) and covariance matrix 
E, i.e. 

f{s,t;z,x) = fo{z- s,x-t), 

where /o is a Gaussian density with mean and covariance matrix E. We assume that the spatial 
dependence is modelled through the matrix S and the temporal dependence is given through 
(T3, leading to the covariance matrix 



af ai2 

ai2 ai I . (2.13) 



S 

o-i/ . „ 



In the following theorem, we calculate a closed form of the bivariate distribution function 
resulting from the setting defined above. The derivation of the bivariate distribution function in 
a purely spatial setting can be found in Padoan, Ribatet and Sisson [30] and we stick closely to 
their notation. The idea of the proof is widely known and for completeness, details are given in 
Appendix El 

Theorem 2.6. With the setting defined above, the max-stable space-time random field 

V{s,t) = \j {^,Uzj-s;xj-t)}, {s,t)eSxT, (2.14) 

has the bivariate distribution function given by 

F{yi,y2) = P{ri{si,ti) < yi,r]{s2,t2) < ^2) 

[ 1 ^ f 2al\og{y2/yi) + ajaikf + uA _ ]_^ / 2ag log(j/i/j/2) + (tJ^W + n' \ ] 

(2.15) 
where h = Si — S2 is the space lag, u = ti — t2 is the time lag and a{h) = {h S~^/i)^'^. 
Note, that if the time lag u equals zero, the formula reduces to 

F{y,,y2) = exp l-lc, f ^ + My^^ _ 1$ f ^ + l°g(^i/^2: 



yi \ 2 a{h) ) y2 \ 2 a{h) 

which is the bivariate distribution of a Gaussian max-stable random field in space as calculated 
in Padoan, Ribatet and Sisson [30]. If the space lag h is zero, the bivariate distribution function 
is given by 

1^/1/, , , . u^W I ^ fl f, , , . u^ 



By comparing the bivariate distributions from the Smith model with those of the approach 
discussed in (j2.6p in Section 12. H we recognize that the functions are the same, if 

6{h, u) = -a(/i)^ + —u^. (2.16) 

In Section HJ where we study a more detailed representation of the function 5, we come back to 
this point. 

3 Pickands dependence function and tail dependence coefficient 

The Pickands dependence function (Pickands [31]) is one measure of tail dependence and is 
related to the so-called exponent measure. A general introduction to exponent measures and 
Pickands dependence function can be found in Beirlant et al. [3]. In particular, the joint distri- 
bution of the max-stable random field can be expressed with the exponent measure V, 

P{r]{si,ti) < yi,r]{s2,t2) < 2/2) = ex.p{-V{yi,y2;6{si - S2,ti - t2))} , 



where in our case V is given through the bivariate distribution function ()2.5p by 

1 / log^ , \ 1 / log2^ , 

yi \2^j6{h,u) J y2 \2y/5{h,u) 

and depends on the space and time lags h and u. 

In the bivariate case, the Pickands dependence function is defined through 

exp {-V{yi, y2, 6{h, u))} = exp <^ - ( \ ] A' 



,yi 2/2/ \yi+y2 

Setting X = yi/{yi + y2), hence 1 - A = 2/2/(2/1 + ^2), we obtain 
A{X; 6ih, u)) = A(l - X)V{X, 1 - A; 6{h, u)) 

A useful summary measure for extremal dependence is the tail-dependence coefficient, which 
goes back to Geffroy [H] and Sibuya [37]. It is defined by 



where F^^ ,n is the generalized inverse of the marginal distribution for fixed location s £ S and 
time point t £ T. For our model this leads to 



X{h,u) = 2{l-^^/ 6{h,u))). (3.1) 

The tail dependence coefficient is a special case of the extremogram introduced in Davis and 
Mikosch [11] (Section 1.4), with the sets A and B defined as (l,oo). The two cases xih,u) = 
and x(^) u) = 1 correspond to the boundary cases of asymptotic independence and complete 
dependence. Thus, if S{h, u) — )■ 0, the marginal components in the bivariate case are completely 
dependent and if S{h, u) — )• 00, the components become independent. In the following section, we 
examine the relationship between the underlying correlation function and the tail dependence 
coefficient. 

4 Possible correlation functions for the underlying space-time 
Gaussian process 

Provided the correlations function of the underlying Gaussian process is sufficiently smooth 
near (0,0), then Assumption 12.11 holds for some sequences Sn and tn- One such condition is 
given below. Throughout this section let h = si — S2 denote the space lag and u = ti — t2 the 
time lag. 

Assumption 4.1. Assume that the correlation function allows for the following expansion 

p{h,u) = l-Ci||/i|ri -Cakr^+Odl/ilpijul^^) 

around (0,0), where < ai,a2 < 2 and Ci,C2 > are constants independent of h and u. 

10 



Remark 4.2. The parameters ai and a2 relate to the smoothness of the sample paths in the 
underlying space-time Gaussian random field. 

Under Assumption 14.11 the scaling sequences in Assumption 12.11 can be chosen as Sn = 
(logn)i/ai and tn = (logn)i/°2. It follows that 

logn{l - p{snh,tnu)) ^ Ci\\h\\°'^ +C2\u\°''^ = 6{h,u), as n — ;> oo. (4.1) 

The condition for the tightness in (|2.3p can be obtained by setting 

Z)((sl,tl),(s2,t2))=max{||sl-S2|r^/^|tl-t2r^/'}, 

which is a metric in M.'^ x [0, oo), since ai,a2 £ (0, 2]. 

Smith's storm profile model can recover a subset of the class of correlation functions specified 
in Assumption 14 . 1 i Choosing ai = 02 = 2 and 

ai2 = 0, aj = al = — , and af = — 

in the Smith model, we find that S{h,u) is of the form given in Assumption 14.11 



6{h, u) = „ o 2^(^2^1 - 2cri2/li/l2 + cr?/ii) + T^U^ 

4(o-f cr| - crfa) 4cr| 



and, hence, has the same finite-dimensional distributions. 

In the following, we analyse several correlation models used in the literature for modelling 
Gaussian random fields in space and time. In recent years, the interest in spatio-temporal correla- 
tion models has been growing significantly; especially in the construction of valid covariance func- 
tions in space and time. A simple way to construct such a model is to take the product of a spatial 
correlation function pi{h) and a temporal correlation function P2{u), i.e. p{h,u) = pi{h)p2{u) 
(see for example Cressie and Huang [U]). Such a model is called separable and Assumption 14.11 
is satisfied, if the spatial and the temporal correlation functions have expansions around zero of 
the form 

pi(/i) = 1 - ciWhr^ + o{\\hr'), P2{u) = 1 - 02^"-' + odur^), 

respectively. 

Example 4.3. A more sophisticated method to obtain covariance models is given on a process- 
based level. An interesting example in this context is presented in Baxevani, Podgorski and 
Rychlik [H [2], who construct spatio-temporal Gaussian random fields in a continuous setup 
using moving averages of spatial random fields over time, given by 



00 
X{s,t)= j f{t-u)^{s,du), 
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where $(-,dti) is a Gaussian random field - valued measure and / is a deterministic kernel 
function. Using the kernel function f{t) = e~^*l|f>o} and the stationary spatial covariance 
model pi{h) = exp | — ||/i|p/C}, one obtains the separable covariance function 

oo 

— oo — oo 

1 1 C llhl|2 

where the temporal dependence is of Ornstein-Uhlenbeck type (see Example 3 in [1]). The 
corresponding correlation function satisfies 

p{h,u) = 1- ^||/if -A|u|+0(||^f |n|). 

Separable space-time models do not allow for any interaction between space and time. Dis- 
advantages of this assumption are pointed out for example in Cressie and Huang [9]. Therefore, 
nonseparable model constructions have been developed. Another approach for combining purely 
spatial and temporal covariance functions leading also to nonseparable covariance models is 
introduced in Ma [261 127] , given in terms of correlation functions by 



oooo 

p{h,u) =11 Pi{hvi)p2{uv2)dG{vi,: 



,V2) 


where G is a bivariate distribution function on [0,oo) x [0,oo). Using the expansions above, it 
follows that 

oo oo 

p{h,u) = l-Ci fv'^'dGi{vi)\\hr^-C2 f v^^dG2{v2)\u\''^ +0{\\h\r\uD, 



where Gi and G2 denote the marginal distributions of G, respectively. From this representation, 
the components in Assumption 14.11 can be defined directly. 

4.1 Gneiting's class of correlation functions 

A more elaborate class of nonseparable, stationary correlation functions is given by Gneiting's 
class [20]. The class of covariance functions is based on completely monotone functions, which 
are defined as functions (p on (0, 00) with existing derivatives of all orders 93'"^ n = 0, 1, . . . and 

(_l)«(^W(t) >o, t>0, n = 0,l,.... 
For our purpose we use a slightly different definition of Gneiting's class. 

Definition 4.4 (Gneiting's class of correlation functions [20]). Let ip : ]R+ — ?• M 6e completely 
monotone and let ip : R+ — )• M 6e a positive function with completely monotone derivative. 
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Further assume that Tp{0) "93(0) = 1, where d is the spatial dimension, and Pi,/32 G (0, 1]. The 
function 

1 f Whf^' \ 

defines a non- separable, isotropic space-time correlation function with p{0,0) = 1. 



Compared to the original definition in [20], we included the parameters /5i and /52, which is 
not a restriction since we can simply change the norms by defining || • H^, and | • 1=,, in terms of 
the old ones through 

||/i||* = ||/i|| , and |u|* = \uf^. 

These new quantities are still norms since /3i,/32 G (0,1]- In the next step, we provide an 
expansion of the correlation function around zero to obtain Assumption 14.11 The following 
proposition generalizes a result by Xue and Xiao [30] (Proposition 6.1). 

Proposition 4.5. Assume that V''(0) 7^ 0. The correlation function taken from the Gneiting 
class satisfies A s sumption \4 -11 with ai = 2/3i, 02 = 2/32 cind 



(4.2) 



(00 00 \ 

I zdF^iz) ^ dF^{z)\ , C2 = ^V(0)"V'(0), 
/ 

where F^p is a non- des creasing bounded function with Fp{0) 7^ and L zdFp{z) < 00. 

Proof. Since the function if is completely monotone, Bernstein's theorem (see for example 
Feller [18], Chapter 13) gives 



00 



^(x) = / e-'^'dF^iz), X > 0. 







From the properties of the correlation function 7/;(0)~ '^93(0) = 1, it follows that V(0) 7^ ^nd 
93(0) 7^ 0. We apply a Taylor expansion to the functions ip{-)~ ''^ and the exponential in the 
representation of ip: 

^{uy^/^ = ^(0)-'^/2 - ^V'(0)-'^/2-V(0)'u + o(u), u^O 

00 00 00 

(p{x) = {1 - xz + o{x))dFp{z) = / dFp{z) - x zdF^{z) + o{x), x —^ 0. 
00 

Using the expansions in the correlation function and replacing u by | u ^ "^ and x by 1 1 h, | p'^^^ / V' ( | if P ) 
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we obtain 



X dF^{z)- zdF^{z)\\hfP^ 



^(0)-'^/2 - --0(O)-'^/2-lyy(o)|^|2/32 ^ o(|u|2/^2) 



+ o{\\hfp^: 



oo 



(oo oo \ 

j zdF^iz)/ I dF^{z)\ \\hf^^ + 0{\\hf^^\u\^^^) 
/ 

= 1 - Ci||/if /^i - C2|n|2/32 + Odl^f /^Mn|2/32), 
where Ci and C2 are defined as in ()4.2p . D 

Remark 4.6. For various choices of ai,a2 € (0,2], i/ie correlation functions in the Gneit- 
ing class have the flexibility to model different levels of smoothness of the underlying Gaussian 
random fields. 

Example 4.7. We illustrate with a specific example, where the functions (p and "0 are taken 
from |20j; namely 

V9(x) = (l + te)"^ tp{x) = {1 + axy , 

where a, 6, i^ > and < 7 < 1. The function y? is the Laplace transform of a gamma probability 
density function with shape z/ > and scale 6 > 0, which has mean bu. The first-order derivative 
of ip at zero is given by ^'(0) = 07. We choose /3i = /32 = 1 leading to ai = 02 = 2 and, thus, a 
mean-square differentiable Gaussian random field. The constants Ci and C2 are given by 

Ci = 6z^, and C2 = —07. 

Figure 14.11 shows contour plots of the correlation function and the resulting tail dependence 
coefficient as in (|3.ip based on different values for a and b with i' = 3/2 and 7 = 1 fixed as a 
function of the space-lag ||^|| and time lag \u\. We see that the tail dependence function exhibits 
virtually the identical pattern of the underlying correlation function under a compression of the 
space-time scale. In particular, the extremal dependence dies out more quickly for large space 
and time lags than for the correlation function. 

In a second step, we simulate random fields in space and time using the above defined 
correlation model with a = b = 0.03, i' = 3/2 and 7 = 1. We start the simulation procedure 
with n = 100 replications of a Gaussian random field with correlation function p{snS, tnu) using 
the simulation routine in the R-package RandomFields by Schlather [3l]. The random fields 
are then transformed to standard Frechet and the pointwise maximum is taken over the 100 
replications. Figure 14.21 shows image and perspective plots (using the R-package fields for 
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visualization) of the simulated random fields for four consecutive time points. Figures 14.31 and 
4.41 show the resulting random fields, if the margins are transformed to standard Gumbel and 
Weibull instead of Prechet distributions. Clearly, in both cases the peaks are not as high as in 
the Frechet case, leading to a smoother appearance of the resulting random field. One still sees 
the isolation of the peaks in the Frechet case, which are well-known from the storm model of 
Smith using a centered Gaussian density for the function /. However, in the other two cases, 
they are not as pronounced. 



Underlying correlation function, a=0.03, b=0.03 



Tali dependence coefficient, a=0.03, b=0.03 



Underlying correlation function, a=0.01 , b=0.05 



Tali dependence coefficient, a=0.01, b=0.05 







10 15 20 



10 15 20 



5 10 15 20 



5 10 15 20 



Under iying correiation function, a =0.0 5, b=0.01 



Tail dependence coefficient, a=0.05, b=0.01 





Figure 4.1: Contour plots for the underlying correlation function (left) and the resulting tail dependence 
coefficient (right) depending on the absolute space lag \\h\\ and time lag |u| for different values of the 
scaling parameters a (time) and b (space), h' = 3/2 and 7=1. 
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Figure 4.2: Simulated max-stable random fields with Frecliet margins (a == 0.03,6 — 0.03, z/ = —3/2, 
7 = 1) for four consecutive time points (from the top to the bottom) with a time lag of one using a grid 
simulation of size 30 x 30. 
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Figure 4.3: Simulated max-stable random fields with Gumbel margins (a = 0.03, b = 0.03, u = —3/2) for 
four eonsecutive time points (from the top to the bottom) with a time lag of one using a grid simulation 
of size 30 X 30. 
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Figure 4.4: Simulated max-stable random fields with Weibull margins (a := 0.03, b = 0.03, u = —3/2) for 
four consecutive time points (from the top to the bottom) with a time lag of one using a grid simulation 
of size 30 X 30. 



18 



4.2 Modelling spatial anisotropy 

The correlation functions of the underlying Gaussian random fields in the previous sections were 
assumed to be spatially isotropic, meaning that the correlation function only depends on the 
absolute space and time lags \\h\\ and \u\. An easy way to introduce spatial anisotropy to a 
model is given by geometric anisotropy, i.e., 

p{h,u) = p{\\Ah\\,\u\), 

where ^ is a transformation matrix. 

In the two dimensional case geometric anisotropy in space can be modelled by a transformation 

matrix A = TR, with rotation and distance matrix where 

I cos a -sina\ /l/omax \ 

I sin a cos a y ' \ ^ l/«min/ 

Geometric anisotropy directly relates to the tail dependence coefficient 



Xih,u) = 2(1 -<^{^5iAh,u))). 

Figure 14.51 compares isotropic and anisotropic correlation functions and the corresponding tail 
dependence coefficients as function of the space lag components h = {hi, /i2)', where the isotropic 
correlation is the same as in Example 14.71 with 7 = 1,1/ = 3/2 and a = b = 0.03. For the 
anisotropic case we choose Omin = l, «max = 3 and a = 45°. It can be seen, that the structure in 
the correlation function translates to the tail dependence coefficient. Corresponding max-stable 
random fields with Frechet margins are shown for four consecutive time points in Figure 14. 6[ 
From the image plots, one clearly sees that the correlation is stronger in one direction. The 
perspective plots show that the isolated peaks are now stretched in one direction. In reality, this 
could correspond to wind speed peaks coming for example from a storm shaped particular in 
this wind direction. 

A more complex way of introducing anisotropy in space is given by the Bernstein class, 
which is introduced in Porcu et al. [32] and revisited in Mateu et al. [29]. The covariance model 
is defined by 

C{h,u) = expl -"^ipidhiDvi - ipti\u\)v2 >dF{vi,V2), 

where F is a bivariate distribution function and ipi,i = 1, . . . ,d and ■0* are positive functions on 
[0, 00) with completely monotone derivatives, also called Bernstein functions. We assume that 
ipi,i = 1, . . . ,d and ipt are standardized, such that tpi{0) = tpt{0) = ^,i = I, ■ ■ ■ ,d. Assumption 
12.11 can directly be derived for the corresponding correlation function. 
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p{h,u) 



C{h,u)/C{0,0) 

d 



i=l 
oo 

l-d 



oo 



l-Y,M\hi\) j I vidF{vi,V2)-MW\) j j V2dF{vi,V2)\ 

*=1 

oooo oooo 

l-d vidF{vi,V2) - I I V2dF{vi,V2)] 

DO 

J oo oo 

l-J2i^-Ci\h,r+o{\hir))jviF„,{vi)-{l-C2\ur+0{\ur))jv2dF,,{v2) 

~^ 

oo oo 

/ vidF^^{vi) - / V2dFt,2{v2)] 


°? d °? d 

1 - Ci y vidF,,{vi) Y, \hir -C2J V2dF,,iv2)\ur + 0{J2 IhiD + 0(1^1°') 



Underlying covariance function, 
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Figure 4.5: Contour plots for covariance functions and tail dependence coefficients depending on the 
space lag components h\ and ft,2 in the isotropic case (top) and for included geometric anisotropy (bottom) 
for a = b = 0.03, amin = 1, Omax = 2 and a = 45°. 
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Figure 4.6: Simulated anisotropic max-stable random fields with Frechet margins from Example (|4.7p 
(a=0.03, b=0.03, u ~ —3/2, 7 = 1) with anisotropic parameters amin = 1, (^max = 3 and a = 45°. 
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5 Conlusion 

The main objective of this paper was to extend concepts of max-stable random fields in space 
to the space-time domain. We extended the idea of constructing max-stable random fields as 
pointwise limits of rescaled and transformed Gaussian random fields, (see Kabluchko et al. [23]) 
In a second step, we extended Smith's storm profile model [38j and calculated the resulting 
bivariate distribution functions. 

We showed that the limit assumption on the correlation function in the underlying Gaussian 
random field relates to the tail dependence coefficient. Several examples of spatio-temporal 
correlation functions and their connection to the tail dependence coefficient have been presented. 
We extended an assumption on the correlation model known from the analysis of extremes of 
stationary Gaussian processes and showed how Gneiting's class of covariance functions [20] 
fits in this context. Visualizations of our results were shown in form of contour plots of the 
underlying correlation functions and the corresponding tail dependence coefficients. In addition, 
we simulated max-stable random fields in space and time using different marginal distributions 
and an anisotropic correlation function. In particular, geometric anisotropy in the underlying 
correlation function lead to directional movements in Smith's storm profile model. 

In a forthcoming paper, composite likelihood and other estimation methods for max-stable 
spatio-temporal random fields are considered. 
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A Derivation of the bivariate distribution function for the space- 
time Smith model 

Proof. (Theorem 12. 6p Since space and time are independent, we can write 

/o(2,x) = /i(z)/2(x), zGM^^xGM 

where /i is the density of a bivariate normal distribution with mean and covariance matrix 
S, and /2 is the denstiy of a normal distribution with mean and variance cj|. Starting from 
equation (I2.12[) with K = 2, we obtain 



oo oo oo 



F{yi,y2) = exp 



Mz,x)\/Mz-h,x-u]^^^^^ 



-oo — oo — oo 
oo oo oo 



exp 



yi / \ y2 

fo{z,x) ( fo{z,x) fo{z-h,x-u) 



-1 



> 



dzdx 



— oo — oo — oo 
oo oo oo 



yi I yi 2/2 

faiz -h,x-u) f fo{z -h,x-u) fo{z, x 



-1 



2/2 



> 



2/2 



2/1 



dzdx 



-oo — oo — oo 



exp{-(I)-(II)} 



a,^//...//fl.{ 



/l(z)^ (h{z)f2{x) ^ f,(z-h)hix-u) 



> ^^-^ '-^^^^ } dzdx 

2/1 2/2 



-oo 

oo 



/2(x)-E 
2/1 



^, MZ)f2{x) ^h{Z-h)Mx-u) 



2/1 



2/2 



dx, 



where Z has a normal density with mean and variance S. Now note that 
lli^lIM > hiZ-h)Mx-u) ^ ^^^^^ ^ ^^^^ _ ^^y, Mx - u) 



2/1 



2/2 



2/2 f2{x) 



^ (2^)-^/2 |S|-i exp |-iz^S-^z| > (2^)-'^/2 \Y.\-^ exp |-^(Z - hfj^-^Z - h)\ 
^ Z^^~'Z < Z^^~'Z - 2Z^^~'h + h^^~'h - 21og (^) - 21og ^^'^"^ " ""^ 



2/1 f2{x-u) 

2/2 f2{x) 



f2{x) 



1 



^ Z'^T.^^h < -h^T,-^h - lo 

" 2 '- V2/2 



2/1 



log 



f2ix - U) 
f2{x) 
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The random variable Z S h =: Z B is normally distributed with mean and variance 



B^^B = h^S^^SS-^/i = h^T^^^h. 



Since /2 is the density of a zero mean normal distribution with variance 0"|, we obtain 



f2 ix-u) \ ^ ]_^ 2 

f2{x) J 2al "' 



logU^^V^ =-__(n2-2nx). 



With a(K) = {h^T.-^hY/'^ and Z'^T.-^h/aih) ~ A/'(0, 1) it follows 

' Z^^^'h a{hl_ log{yi/y,) 1 , _ \ 

a{h) - 2 a{h) ^ 2ala{hy ') 



' " a(/i) 2(j|a(/i) y ■ 



Altogether, for independent random variables N and X with A^ standard normally distributed 
and X normally distributed with mean and variance a^, it holds 

oo 

yi J V 2 a(/i) 2a^a{h) cr^a[h) J 

—oo 

yi V a3a{h)a3~ 2 a(/i) 2ala{h) ) 

/ a(h) 1^ log(3/2/j/i) _i_ it^ \ / \ 

]_^ / ^~ + a{h) + 2^^^(7r) \ ^ J_^ / 2a| log(j/2M) + ^lajhf + u^ \ 



since N + u/{a3a{h)){X/a3) is normally distributed with mean and variance l + v?' /{a'^afji)^). 
Analogously to (I), using the substitution Z — t- Z + h, we obtain the second term (II) in 
dUSD. D 
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